realizationsCount = 0;
for k = 1:1:100
	
%parameter: realizations 
realizations = 1000;

%parameters Xa, Xb: X min and max values
Xa = -0.5;
Xb = 0.5;

%parameters Ya, Yb: Y min and max values
Ya = -0.5;
Yb = 0.5;

%constant: Int. Represents the analytical exact value of the calculation
Int = 0.148496;

YBorderTerm = 1 / (Yb - Ya);
XBorderTerm = 1 / (Xb - Xa);

ihat = 0;
keep = 1;
i =1;
while keep == 1
	
	x = rand() - 0.5;
	y = rand() - 0.5;
	
	squareSum = (x^2 + y^2);

	f = squareSum * e ^ (- squareSum / 2) / (YBorderTerm * XBorderTerm);
	ihat = ihat + f;
	
	value(i) = ihat / i;
	
%

media = 0;
for p = 1:1:i
	media = media +  value(p);
end

media = media / i;
suma = 0;
j = 2;
tenpercent = 0;
desv = 1000;
while j <= i
	valor = (value(j) - media)^2;
	suma = suma + valor;

	desv = sqrt(suma / (j - 1));
	
	j = j + 1;
end

%	if abs(value(j-1) - Int)  < 0.1 * Int 
	if desv < 0.01 * Int
		keep = 0;
		realizationsCount = realizationsCount + i + 1;	
		desv
		N = i + 1
		value(j -1)
	end
%	
i = i + 1;
end

end
m = realizationsCount / 100

